/*
 * System.cpp
 *
 *  Created on: 22 Aug 2011
 *      Author: Allan
 */

#include "System.h"

// C++ includes
#include <algorithm>

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "Utilities/Algorithms.h"
#include "Utilities/PhaseUtils.h"
#include "Utilities/ReactionUtils.h"
#include "Utilities/SpeciesUtils.h"

System::System()
{}

System::System(const vector<string>& species)
{
	vector<string> systemSpecies = Union(DetermineAqueousSpecies(species), species);
	
	vector<string> eSpecies = FilterSecondarySpecies(systemSpecies);
	
	multiphase = CreateMultiphase(systemSpecies);
	eReactions = CreateReactionSystem(eSpecies, multiphase.GetSpeciesMap());
	
	InitialiazeReactionMap();
}

System::System(const vector<string>& species, const vector<string>& kSpecies)
{
	vector<string> systemSpecies = Union(Union(DetermineAqueousSpecies(species), species), kSpecies);
	
	vector<string> eSpecies = Difference(FilterSecondarySpecies(systemSpecies), kSpecies);
	
	multiphase = CreateMultiphase(systemSpecies);
	eReactions = CreateReactionSystem(eSpecies, multiphase.GetSpeciesMap());
	kReactions = CreateReactionSystem(kSpecies, multiphase.GetSpeciesMap());
	
	InitialiazeReactionMap();
}

System::System(const Multiphase& multiphase, const ReactionSystem& eReactions) :
multiphase(multiphase), eReactions(eReactions)
{
	InitialiazeReactionMap();
}
	
System::System(const Multiphase& multiphase, const ReactionSystem& eReactions, const ReactionSystem& kReactions) :
multiphase(multiphase), eReactions(eReactions), kReactions(kReactions)
{
	InitialiazeReactionMap();
}
	
const Multiphase& System::GetMultiphase() const
{
	return multiphase;
}

const ReactionSystem& System::GetEquilibriumReactions() const
{
	return eReactions;
}

const ReactionSystem& System::GetKineticReactions() const
{
	return kReactions;
}

const Phase& System::GetPhase(const string& phasename) const
{
	auto checkName = [&](const Phase& phase) { return phase.GetName() == phasename; };
	
	const vector<Phase>& phases = multiphase.GetPhases();
	
	return *find_if(phases.begin(), phases.end(), checkName);
}

const VectorXd System::GetPhaseComposition(const string& phasename, const VectorXd& n) const
{
	unsigned start = 0;
	unsigned size = 0;
	
	BOOST_FOREACH(const Phase& phase, multiphase.GetPhases())
	{
		if(phase.GetName() == phasename)
			size = phase.GetSpecies().size();
		else
			start += phase.GetSpecies().size();
	}
	
	return n.segment(start, size); 
}

const vector<string> System::GetPrimarySpecies() const
{
	return Union(eReactions.Products(), kReactions.Products());
}

const vector<string> System::GetEquilibriumSpecies() const
{
	return eReactions.Reactants();
}

const vector<string> System::GetKineticSpecies() const
{
	return kReactions.Reactants();
}

const Index	System::operator[](const string& species) const
{
	return multiphase.GetSpeciesMap()[species];
}
	
const vector<Index> System::operator[](const vector<string>& species) const
{
	return multiphase.GetSpeciesMap()[species];
}
	
void System::Print() const
{
	cout << "Multiphase System: " << endl;
	cout << multiphase << endl << endl;
	
	cout << "Equilibrium-Controlled Reactions: " << endl;
	cout << eReactions << endl << endl;
	
	cout << "Kinetically-Controlled Reactions: " << endl;
	cout << kReactions << endl << endl;
}

void System::InitialiazeReactionMap()
{
	BOOST_FOREACH(const string& species, GetPrimarySpecies())
		reactionMap[species] = CreateBasisReaction(species, multiphase.GetSpeciesMap());
	
	BOOST_FOREACH(const Reaction& reaction, eReactions.GetReactions())
		reactionMap[reaction.GetReactant().name] = reaction;
	
	BOOST_FOREACH(const Reaction& reaction, kReactions.GetReactions())
		reactionMap[reaction.GetReactant().name] = reaction;
}
